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Abstract 

This  project  focused  on  DDDAS -motivated  developments  in  support  of  space  weather  monitoring  and 
prediction.  The  project  involved  four  interrelated  tasks  relating  to  physics-driven  adaptive  modeling,  adap¬ 
tive  data  assimilation  with  input  reconstruction,  event-based  sensor  reconfiguration,  and  optimization  of 
scheduling.  For  data  assimilation,  the  emphasis  has  been  on  model  refinement.  The  problem  of  estimat¬ 
ing  the  eddy  diffusion  coefficient  using  total  electron  content  measurements  has  led  to  new  techniques  for 
determining  the  essential  modeling  details  needed  by  the  retrospective  cost  model  refinement  technique. 
For  spacecraft  design,  multidisciplinary  optimization  design  techniques  were  applied  to  the  design  of  small 
satellites  accounting  for  multiple  vehicle  subsystems.  For  download  scheduling,  optimization  techniques 
were  used  to  account  for  multiple  spacecraft  and  ground  stations. 
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1  Project  Objectives 


Space  weather  can  have  a  catastrophic  effect  on  operational  satellites  and  ground  systems.  The  ability  to 
monitor  space  weather  is  thus  important  to  both  the  U.S.  civilian  and  military  infrastructure.  With  this 
motivation  in  mind,  this  project  addresses  aspects  of  each  of  the  four  DDDAS  science  and  technology 
Frontiers.  For  the  Applications  Modeling  Frontier,  we  consider  physics-driven  adaptive  modeling ,  where 
current  conditions  elicit  the  need  for  an  appropriate  physics  model.  For  the  Mathematical  and  Statistical 
Algorithms  Frontier,  we  consider  adaptive  data  assimilation  with  input  reconstruction,  where  the  ability 
to  reconstruct  inputs  can  enhance  the  ability  to  estimate  the  state  of  a  system  and  improve  the  accuracy  of 
the  underlying  model.  For  the  Application  Measurement  Systems  and  Methods  Frontier  we  consider  event- 
based  sensor  reconfiguration,  where  a  measurement  system  modifies  its  functionality  based  on  operational 
priorities.  For  the  Advances  in  Systems  Software  Frontier  we  consider  hierarchical  closed-loop  scheduling, 
where  communication  links  and  their  constraints  arc  accounted  for  in  both  autonomous  agent-based  and 
centralized  software  systems. 


2  Summary  of  Main  Results 

2.1  Model  Refinement  and  Input  Reconstruction 

For  modeling  the  ionosphere-troposphere,  we  use  the  Global  Ionosphere-Thermosphere  Model  [77],  which 
is  a  CFD  code  with  atmospheric  chemistry,  as  the  basis  for  data  assimilation,  input  estimation,  and  subsystem 
identification.  These  objectives  were  addressed  by  applying  a  specialized  technique  that  is  being  developed 
under  this  project.  This  technique  is  called  retrospective-cost-based  adaptive  input  and  state  estimation 
(RCAISE).  We  applied  RCAISE  to  GITM  in  order  to  estimate  unknown  drivers  (inputs)  to  the  ionosphere- 
thermosphere  [2,24].  This  approach  complements  ensemble-based  methods  as  in  [66],  which  applies  the 
ensemble  code  DART  to  GITM,  since  the  main  goal  of  RCAISE  is  to  estimate  the  unknown  input  rather  than 
obtain  estimates  of  all  model  states.  Since  RCAISE  is  not  an  ensemble  code,  it  is  highly  computationally 
efficient  compared  to  DART,  but,  unlike  DART,  it  provides  best  estimates  rather  than  a  probability  density 
function.  A  related  approach  was  previously  demonstrated  in  [25]  to  identify  a  subsystem  cooling  model 
whose  dynamics  are  inaccessible  in  the  sense  that  neither  the  input  nor  the  output  of  the  cooling  subsystem 
dynamics  arc  measured. 

Much  of  this  effort  has  focused  on  determining  the  essential  modeling  details  required  by  RCAISE  for 
both  input  estimation  and  model  refinement.  For  linear  systems,  this  information  is  well  understood.  How¬ 
ever,  for  application  to  space  weather  modeling,  GITM  is  both  extremely  nonlinear  and  high  dimensional. 
Our  goal  is  thus  to  extract  the  essential  modeling  details  from  numerical  simulations.  For  earlier  studies, 
numerical  testing  was  used  to  determine  this  information.  However,  an  additional  goal  was  to  estimate 
the  eddy  diffusion  coefficient  (EDC)  using  total  electron  content  (TEC)  data.  TEC  data  is  available  from 
ground  stations  located  worldwide,  and  thus  RCAISE  is  able  to  work  with  multiple  measurements.  Numeri¬ 
cal  testing  has  shown  that  estimating  this  modeling  information  is  computationally  expensive,  and  thus  more 
efficient  methods  arc  needed. 

In  previous  work,  RCMR  was  used  in  [25]  to  estimate  NOx  cooling  profile  using  GITM  and  simulated 
spacecraft  measurements.  This  study  was  limited  to  a  one-dimensional  version  of  GITM.  RCMR  was  sub¬ 
sequently  used  in  [18]  to  estimate  the  photoelectron  heating  coefficient  in  a  fully  three-dimensional  version 
of  GITM  using  both  simulated  and  real  satellite  measurements.  The  applications  in  [25]  and  [18]  demon¬ 
strated  the  ability  of  RCMR  to  adaptively  refine  a  nominal  model  by  iteratively  updating  an  estimate  of  an 
unknown  parameter  modeled  by  an  adaptive  subsystem  model.  This  parameter  is  modeled  as  an  unknown 
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Figure  1:  Retrospective  Cost  Model  Refinement  (RCMR)  architecture.  The  main  physical  system  is  inter¬ 
connected  by  means  of  feedback  with  an  unknown  subsystem  that  models  the  uncertainty  in  the  physical 

system.  RCMR  uses  the  model  output  error  z  =  yo  —  yo  to  adjust  the  adaptive  subsystem  model  to  estimate 
the  unknown  subsystem.  The  update  algorithm  is  based  on  retrospective  cost  optimization. 


subsystem,  as  shown  in  Figure  1 ,  where  the  input  signal  y  and  output  signal  u  of  the  unknown  subsystem 
arc  not  measured  and  thus  cannot  be  used  to  estimate  the  unknown  subsystem. 

Alternatively,  for  the  objective  of  estimating  an  unknown  external  input,  RCAISE  was  used  in  [2]  to 
estimate  the  unknown  driver  F10.7  in  a  fully  three-dimensional  version  of  GITM  using  both  simulated  and 
real  satellite  measurements.  In  this  case,  as  shown  in  Figure  2,  the  adaptive  driver  estimator  is  updated  to 
obtain  an  estimate  w  of  the  unknown  driver  w.  These  studies  demonstrated  that  RCMR  and  RCAISE  can 
effectively  use  data  parameter,  state,  and  input  estimation  within  the  context  of  a  highly  nonlinear,  large- 
scale  model  consisting  of  thousands  of  state  variables. 

2.2  Estimation  of  the  Eddy  Diffusion  Coefficient 

In  this  section  we  briefly  describe  work  where  the  goal  is  to  estimate  the  eddy  diffusion  coefficient  (EDC) 
of  the  thermosphere.  The  drag  force  felt  by  low-Earth  orbiting  objects  is  linearly  proportional  to  the  mass 
density  of  the  thermosphere.  Uncertainties  in  thermospheric  mass  density  variations  arc  the  major  limit¬ 
ing  factor  for  precise  low-Earth  orbit  determination.  The  perturbation  of  the  thermospheric  mass  density  is 
strongly  controlled  by  the  energy  deposited  into  the  upper  atmosphere.  The  difference  in  the  thermosphere 
mass  density  responses  to  different  sources  of  energy  input  has  been  investigated,  and  the  spatial  and  tem¬ 
poral  variations  of  the  thermospheric  mass  density  during  a  series  of  idealized  substorms  were  studied  using 
the  Global  Ionosphere  Thermosphere  Model  (GITM)  [64],  The  mass  density  response  to  different  types  of 
energy  inputs  was  shown  to  have  strong  local  time  dependence.  In  addition,  the  thermosphere  mass  density 
response  to  different  sources  of  energy  input  is  a  slightly  nonlinearly  system,  and  the  non-linearity  grows 
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Figure  2:  Retrospective  Cost  Adaptive  Input  and  State  Estimation  (RCAISE)  architecture.  RCAISE  uses  the 
estimation  error  z  =  y o  —  yo  to  adjust  the  adaptive  driver  estimator,  which  is  used  to  construct  an  estimate 
w  of  the  unknown  driver  w.  The  update  algorithm  is  based  on  retrospective  cost  optimization. 


with  the  level  of  energy  input. 

The  Eddy  diffusion  coefficient  represents  the  mixing  of  the  lower  thermosphere  composition.  By  control¬ 
ling  the  thermo-conductions  and  the  altitudes  where  the  species  begin  to  diffuse  according  to  their  particular 
scale  height,  variations  in  the  Eddy  diffusion  coefficient  strongly  affect  the  upper  thermosphere,  and  hence 
the  ionosphere  total  electron  content  (TEC),  which  is  a  key  parameter  in  the  ionosphere.  A  goal  of  this  work 
is  to  understand  how  EDC  affects  TEC.  The  effect  of  EDC  variations  on  TEC  was  investigated  in  detail 
under  various  solar  conditions,  including  the  dependence  of  EDC  on  the  F10.7  index. 

The  thermosphere  composition  change  has  a  strong  effect  on  the  thermosphere  mass  density  response, 
especially  in  the  oxygen/helium  transition  region;  however,  most  physical  models  do  not  include  the  effect 
of  helium.  To  better  represent  the  thermosphere  mass  density  perturbations  during  geomagnetic  storms,  the 
helium  component  is  implemented  in  the  GITM  model.  In  ongoing  research,  the  morphology  of  helium 
simulated  by  the  GITM  model  is  being  investigated  and  compared  with  measurement  from  the  Atmosphere 
Explorer  (AE)  satellite. 

The  parameter  EDC  plays  a  key  role  in  GITM  and  many  CFD  (computational  fluid  dynamics)  codes. 
In  particular,  EDC  in  GITM  models  the  turbulent  mixing  in  the  upper  atmosphere.  According  to  the  mass 
continuity  and  the  momentum  equations,  the  altitude  profile  of  the  neutral  constituents  changes  from  fully 
mixing  at  lower  altitudes,  where  turbulent  mixing  prevails,  to  molecular  diffusion  at  higher  altitudes  (above 
110  km).  The  value  of  EDC  represents  the  intensity  of  the  turbulent  mixing,  which  is  a  key  factor  in 
controlling  the  upper  atmosphere  composition  profiles. 

In  practice,  it  is  challenging  to  estimate  EDC  from  ground-based  and  satellite  measurements  [39,58, 
83],  and  efforts  to  estimate  the  value  of  EDC  based  on  theoretical  analysis  have  been  made  by  [1,41,73] 
Estimates  of  the  EDC  value  based  on  measurements  and/or  theoretical  modeling  range  from  25  to  1000 
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m2/sec.  In  view  of  these  diverse  estimates,  it  is  desirable  to  develop  techniques  that  can  estimate  EDC 
based  on  the  available  data. 


Figure  3:  Simulated  Ground  Locations  for  Total  Electron  Count  Measurements.  Each  ”+”  denotes  the 
location  of  a  simulated  ground  location  for  measurements  of  total  electron  count  (TEC).  Simulated  data 
from  these  locations  is  used  by  RCMR  to  estimate  the  eddy  diffusion  coefficient  in  the  thermosphere. 


RCMR  and  RCAISE  depend  on  modeling  information  about  the  system  being  considered;  this  modeling 
information  can  be  viewed  as  a  reduced-order  model,  or  as  tuning  parameters,  since  often  only  a  small 
number  of  parameters  are  needed.  For  a  low-dimensional  linear  system,  the  required  modeling  information 
can  be  extracted  in  the  form  of  the  impulse  response.  For  DDDAS  applications,  however,  our  interest  is  in 
high-dimensional  nonlinear  systems,  such  as  the  thermosphere  as  modeled  by  GITM,  where  the  number  of 
states  may  be  as  large  as  10'.  GITM  is  implemented  by  a  FORTRAN  code  that  captures  diverse  physics 
including  fluid  dynamics,  thermodynamics,  chemical  kinetics,  and  electrodynamics.  For  this  application,  no 
analytical  model  is  available,  and  it  is  not  possible  to  analytically  extract  the  required  modeling  information 
due  to  the  complexity  of  the  physics.  Therefore,  we  typically  implement  RCMR  and  RCAISE  using  a 
small  number  of  tuning  parameters  determined  by  numerical  testing.  As  we  have  found  in  relation  to  the 
estimation  of  EDC  using  TEC  measurements,  this  approach  is  inefficient. 

What  is  needed  to  make  RCMR/RCAISE  a  truly  practical  tool  is  a  reliable,  systematic,  and  easily  imple- 
mentable  technique  for  obtaining  the  tuning  parameters  required  by  RCMR/RCAISE.  One  approach,  which 
we  have  tested  on  low-order  examples,  involves  iterative  refinement  of  the  tuning  parameters  that  comprise 
the  filter  Gf.  For  a  fixed  data  window,  initial  tuning  parameters  are  chosen  and  used  to  obtain  a  subsys¬ 
tem  model.  This  subsystem  model  is  merged  with  main  system  model,  and  the  resulting  “closed-loop” 
impulse  response  is  computed.  The  impulse  response  parameters  are  then  used  as  tuning  parameters  within 
an  updated  filter  Gf  for  re-estimating  the  unknown  subsystem,  and  the  process  repeats  using  the  same  data 
window  used  in  previous  iterations.  Tests  have  shown  that  this  iterative  filter  refinement  approach  provides 
fast  response  from  a  possibly  poor  initial  choice  of  Gf  and  is  efficient  in  the  sense  that  only  limited  data 
window  is  needed. 

To  illustrate  this  technique,  consider  a  2nd  order  linear  system  where  An  is  uncertain.  We  assume  that 
the  nominal  model  of  A  is  such  that  An  —  An  =  AAn  =  1.  We  construct  Gf  from  the  impulse  response  of 
transfer  function  from  u  to  yo-  In  subsequent  iterations,  Gf  is  refined  using  the  updated  value  of  A.  Figure 
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(a) 


Figure  4:  (a)  RCMR  estimate  of  A  An  with  iterative  filter  refinement.  The  filter  Gf  is  refined  after  each 
iteration  based  on  the  updated  estimate  of  AAn.  Note  that  the  data  window  for  each  iteration  consists  of 
only  10  time  steps,  RCMR  is  switched  on  after  5  time  steps,  and  one-step  estimation  of  the  estimate  of  AAn 
is  achieved  at  the  3rd  iteration,  (b)  Estimate  of  AAn  at  the  end  of  each  iteration.  After  each  iteration,  Gf  is 
refined  using  the  updated  estimate,  and  is  used  in  the  next  iteration. 


4(a)  shows  the  estimate  of  AA\ \  in  each  iteration.  Figure  4(b)  shows  the  estimate  of  A4n  at  the  end  of 
each  iteration. 

This  Gf  refinement  technique  has  the  added  advantage  that  refinement  of  Gf  is  implemented  on  the  same 
data  set,  and  thus  only  a  small  amount  of  data  is  required  to  increase  the  rate  of  convergence.  The  final  refined 
filter  Gf  is  then  used  on  the  remaining  data  set.  For  nonlinear  systems,  this  technique  can  simplify  the  numer¬ 
ical  testing  procedure  for  constructing  Gf.  If  successful  on  large-scale  applications  such  as  GITM,  this  iter¬ 
ative  filter  refinement  technique  will  significantly  simplify  the  implementation  of  RCMR/RCAISE  on  new 
and  challenging  applications.  In  particular,  this  technique  will  facilitate  the  application  of  RCMR/RCAISE 
with  greater  assurance  of  success. 

2.3  Cubesat  Technology 

An  additional  component  of  this  project  was  the  development  of  post-launch  calibration  techniques  for 
spacecraft  that  collect  data  for  atmospheric  modeling.  This  effort  ties  into  the  development  of  Cubesats 
at  the  University  of  Michigan,  several  of  which  are  currently  in  orbit  [11,22,78].  New  techniques  were 
developed  for  calibrating  on-board  sensors  after  the  spacecraft  reaches  its  specified  orbit.  In  particular, 
low-cost  photoelectric  and  magnetic-field  sensors  typically  lose  calibration  due  to  the  launch  and  space 
environment.  Through  a  combination  of  physics  modeling  and  data  analysis,  it  is  shown  in  [89]  that  the 
accuracy  of  on-board  sensors  can  be  enhanced  through  ground-based  re-calibration.  This  re-calibration 
accounts  for  time-varying  electrical  and  magnetic  fields,  the  Earth’s  spatially  varying  magnetic  field,  and 
sensor  degradation. 

Additional  research  has  focused  on  design  studies  for  spacecraft  configurations.  In  particular,  multidis¬ 
ciplinary  design  optimization  techniques  were  applied  in  [5 1]  to  optimize  multiple  spacecraft  subsystems, 
including  power  and  communications  subsystems.  In  addition,  a  study  of  photovoltaic  power  generation 
constraints  due  to  spacecraft  solar  panel  geometry,  orbital  effects,  and  self-induced  shadowing  is  described 
in  [61]. 
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2.4  Scheduling  for  Satellite  Data  Downloading 


The  last  component  of  this  project  addresses  the  challenges  that  arise  from  transferring  data  from  a  con¬ 
stellation  of  satellites.  Data  transfer  is  limited  due  to  power  and  energy  constraints,  orbit  trajectory,  and 
the  location,  gain,  and  field  of  view  of  available  ground-based  stations.  This  leads  to  a  capacity-constrained 
scheduling  problem.  These  issues  were  addressed  by  developing  real-time  scheduling  algorithms  for  model¬ 
ing  and  optimizing  space  networks  to  optimize  communication  capacity  [86,  88].  Research  includes  detailed 
investigations  of  download-scheduling  techniques  for  configurations  involving  multiple  satellites  and  ground 
stations  [19,62,  87]. 

2.5  Student  Training 

This  project  partially  supported  several  students  who  received  the  Ph.D.  degree,  namely,  Sara  Spangelo, 
John  Springmann,  and  Asad  Ali,  as  well  as  postdoc  Angeline  Burrell.  Continuing  students  include  Ankit 
Goel  in  Aerospace  Engineering  as  well  as  postdoc  XianJing  Liu.  Brian  Lemay  and  Jeremy  Castaing  arc 
expected  to  complete  their  Ph.D.  degrees  in  2016  and  2017,  respectively. 

Two  students  supported  by  this  project  were  recognized  for  their  technical  contributions.  John  Spring¬ 
mann  was  the  winner  of  the  Student  Paper  Competition  at  the  AIAA  2013  Small  Satellite  Conference 
for  [90].  Jeremy  Castaing  received  honorable  mention  in  the  Student  Paper  Competition  at  the  AIAA  2014 
Small  Satellite  Conference  for  [19]. 


3  Background:  Space- Weather  Monitoring  for  Safety  and  Reliability 

The  near-Earth  environment,  in  particular-,  the  ionosphere-thermosphere,  changes  in  response  to  solar  phe¬ 
nomena.  These  phenomena  include  massive  ejections  of  charged  particles,  which  can  disrupt  the  upper 
atmosphere  and  magnetic  field  of  the  Earth.  These  changes  can  have  catastrophic  effects  of  various  types. 
First,  changes  in  electrical  properties  affect  the  transmission  of  radio  waves,  which  degrades  the  accuracy 
of  GPS.  Next,  changes  in  density  and  wind  speed  affect  the  orbits  of  both  operational  satellites  and  de¬ 
bris.  These  orbital  changes  make  it  difficult  to  predict  the  possibility  of  collisions  and  take  evasive  action. 
Collisions  have  destructive  consequences,  and  lead  to  more  debris,  which  further  exacerbates  the  problem. 
Finally,  space  weather  can  induce  ground  currents,  which  can  destroy  electrical  transformers.  The  arti¬ 
cle  [53]  describes  the  potentially  catastrophic  consequences  of  such  events,  which  have  occurred  several 
times  during  the  last  century  and  whose  reoccurrence  is  thought  to  be  inevitable. 

The  accuracy  of  predictions  of  satellite  trajectories  through  the  upper  atmosphere  depends  on  the  ac¬ 
curacy  of  the  drag  model.  The  drag  force  F,/  is  modeled  by  F,/  =  —  \ACdp\\Vs  —  ||2,  where  A  is  the 

projected  surface  area  of  the  satellite,  C,/  is  the  drag  coefficient,  p  is  the  density,  Vs  is  the  velocity  of  the 
satellite,  and  Vw  is  the  velocity  of  the  upper  atmospheric  wind  that  the  satellite  encounters.  In  the  upper 
atmosphere,  the  wind  is  typically  200-400  m/s,  while  the  satellite  orbital  velocity  is  on  the  order  of  7500 
m/s;  therefore,  the  wind  is  typically  ignored.  During  strong  upper  atmospheric  storms,  however,  the  winds 
can  increase  beyond  1000  m/s.  The  density  p  is  even  more  variable.  For  example,  from  night  to  day,  p 
can  increase  by  a  factor  of  two.  This  means  that  as  the  satellite  orbits  (in  90  minutes),  it  traverses  from 
the  sunlit  dayside  to  the  colder  nightside,  encountering  very  different  drag  conditions.  Furthermore,  from 
solar  minimum  conditions,  when  the  sun  is  dimmest  in  X-ray  wavelengths,  to  solar  maximum  (5-6  years 
later),  when  it  is  brightest  in  those  wavelengths,  the  nominal  density  can  increase  by  more  than  an  order  of 
magnitude.  Density  can  also  increase  dramatically  on  short  time-scales;  for  example,  during  a  storm  p  has 
significant  localized  structure  and,  globally  averaged,  can  increase  by  a  factor  of  two  or  three  in  just  a  few 
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hours  (3-4  orbital  periods).  These  types  of  storms  cause  the  Joint  Space  Operations  Center  to  miscalculate 
the  trajectories  of  thousands  of  objects,  and  days  of  effort  arc  needed  to  re-acquire  their  orbital  positions. 

Fig.  5  illustrates  the  variation  in  p  observed  by  the  CHAMP  satellite  [10]  over  four  days,  during  which 
two  large  storms  occurred.  The  upper  plot  shows  1 -minute  averaged  data,  while  the  lower  plot  shows  data 
averaged  over  a  90-minute  orbit.  The  maximum  density  observed  by  the  satellite  is  roughly  28  x  10“ 12 
kg/m3,  while  the  smallest  density  is  about  4  x  10“ 12  kg/m3,  that  is,  a  factor  of  seven  smaller.  The  variability 
observed  over  the  orbit  in  the  upper  plot  is  a  result  of  the  satellite  passing  through  the  warm  summer  dayside 
into  the  cold  winter  nightside.  When  the  data  arc  averaged  over  an  orbit,  increases  during  the  storm  times 
arc  apparent.  The  rises  arc  very  fast,  and  p  approximately  doubles  over  the  course  of  a  few  hours. 

Empirical  models  of  the  upper  atmospheric  density  arc  used  to  predict  orbital  trajectories.  Since  these 
empirical  models  were  found  to  lack  the  ability  to  accurately  estimate  p,  the  Department  of  Defense  launched 
approximately  200  spheres,  which  simply  orbit  the  Earth  and  arc  affected  by  the  drag.  Their  position  is 
determined  by  radar  twice  a  day,  so  that  the  empirical  models  can  be  “corrected”  by  scaling  to  compensate 
for  the  global  error.  Although  these  empirical  models  arc  easy  to  run,  they  arc  crude  and  have  limited 
accuracy. 

The  most  proactive  approach  to  dealing  with  the  effects  of  space  weather  is  to  sense  disturbances  at  the 
earliest  possible  stage  and  take  evasive  action.  Long-term  (up  to  three-day)  prediction  requires  measure¬ 
ments  of  solar  activity  collected  by  spacecraft  in  solar  orbit  as  well  as  spacecraft  at  the  Lagrangian  point  that 
monitor  the  solar  wind.  These  sensors  provide  advance  knowledge  of  solar  storms,  which  can  then  be  used 
to  predict  the  impact  on  the  ionosphere-thermosphere.  The  project  that  we  propose  focuses  on  one  com¬ 
ponent  of  the  space  weather  problem,  namely,  monitoring  of  the  ionosphere-thermosphere.  As  we  describe 
below,  this  objective  requires  multiple  transformative  advances  in  DDDAS. 


4  Relevance  of  DDDAS  to  Space  Weather  Monitoring 

The  long-term  motivation  for  this  DDDAS  project  is  to  develop  a  constellation  of  Cubesats  that  can  provide 
data  to  dynamically  monitor  the  ionosphere-thermosphere.  A  constellation  of  approximately  40  Cubesats  in 
different  orbits  would  provide  the  ability  to  monitor  the  atmosphere  in  a  global  24/7  manner.  Each  Cube- 
sat  would  be  equipped  with  a  suite  of  sensors  for  monitoring  properties  of  the  ionosphere-thermosphere. 
Data  would  be  stored  on-board  and  communicated  to  the  ground  where  it  could  be  used  for  modeling  and 
prediction. 

The  operation  of  this  constellation,  however,  requires  multiple  transformative  advances  in  DDDAS.  For 
starters,  each  Cubesat  has  limited  energy  and  power,  which  constrains  the  ability  of  each  Cubesat  to  send 
and  receive  data.  The  power  needed  for  communication  depends  partially  on  the  ground  station,  and  ground 
stations  vary  widely  in  their  ability  to  send  and  receive  data  as  well  as  their  availability,  which  can  change 
from  hour  to  hour.  Further  complicating  the  communication  problem  is  the  fact  that  each  Cubesat  is  in  a 
different  orbit  and  has  only  partial  information  about  the  availability  of  ground  stations  along  its  path.  At 
the  same  time,  the  ground  controller  has  limited  ability  to  communicate  with  each  Cubesat  to  update  it  on 
ground-station  availability. 

While  these  issues  suggest  the  dynamic  nature  of  the  problem  of  the  data  collection  problem,  the  DDDAS 
nature  of  the  problem  becomes  much  more  complex  when  the  dynamics  of  the  atmosphere  arc  taken  into 
account.  Consider,  for  example,  the  following  scenario.  Each  Cubesat  has  an  on-hoard  GPS  system,  and  it 
also  has  some  computational  ability  that  it  can  use  to  predict  its  own  location.  When  a  Cubesat’s  predicted 
location  varies  significantly  from  its  measured  position,  it  can  infer  the  presence  of  perturbations  to  the 
atmospheric  density.  The  Cubesat  can  then  begin  to  collect  data  at  higher  resolution  than  normal  and,  at  the 
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same  time,  must  determine  an  efficient  means  for  transmitting  its  data  to  a  ground  station.  The  transmitted 
data  arc  ultimately  collected  centrally,  where  it  can  be  used  for  data  assimilation.  The  results  of  the  data 
assimilation,  which  must  be  performed  at  nearly  a  real-time  rate,  can  be  used  to  determine  whether  the 
disruption  to  the  ionosphere-thermosphere  is  significant.  If  so,  additional  Cubesats  can  be  enlisted  to  begin 
data  collection  in  an  optimal  configuration  and  at  an  enhanced  rate,  while  the  downlinking  for  each  Cubesat 
must  be  scheduled  “on  the  fly”  based  on  priorities  set  by  the  data  assimilation. 

This  is  only  one  scenario  demonstrating  how  sensing,  scheduling,  data  assimilation,  and  analysis  of  the 
results  must  work  together  to  monitor  the  ionosphere.  The  key  point  is  that  each  component  of  this  system 
is  dynamic  and  data  driven,  and  that  data  assimilation  must  interact  with  data  collection  and  vice  versa  in 
real  time.  This  application  thus  reflects  the  goals,  challenges,  and  opportunities  of  the  DDDAS  paradigm. 

DDDAS  seeks  transformative  advances  in  four  frontiers,  namely.  Applications  Modeling,  Mathemati¬ 
cal  and  Statistical  Algorithms,  Application  Measurement  Systems  and  Methods,  and  Advances  in  Systems 
Software.  Research  in  each  frontier  is  motivated  by  the  dynamic  data-driven  aspects  of  the  Cubesat  con¬ 
stellation,  data  assimilation,  and  space  weather  conditions,  all  of  which  are  essential  to  the  overarching  goal 
of  monitoring  space  weather  in  real  time.  To  achieve  the  required  advances  in  each  frontier,  the  following 
sections  describe  relevant  research  in  each  area.  Each  section  describes  the  motivation  for  DDDAS-related 
research  within  the  context  of  space  weather  monitoring  along  with  recent  research  advances. 


5  DDDAS  Frontier:  Applications  Modeling 

5.1  Context  and  Motivation  within  DDDAS 

The  foundation  for  space  weather  monitoring  is  a  first  principles  model  of  the  upper  atmosphere,  namely, 
the  Global  Ionosphere-Thermosphere  Model  (GITM)  [77],  which  provides  the  basis  for  data  assimilation 
algorithms.  GITM  is  a  3-dimensional  spherical  code  that  solves  the  Navier-Stokes  equations  for  the  thermo¬ 
sphere.  These  types  of  models  work  better  than  empirical  (ad  hoc  “correction”)  models  because  they  capture 
the  dynamics  of  the  system  instead  of  snapshots  of  steady-state  solutions,  which  arc  what  most  empirical 
models  provide.  Furthermore,  first  principles  models,  such  as  GITM,  model  the  winds  that  can  influence  the 
drag.  In  order  to  accurately  predict  p  and  Vw  in  the  upper  atmosphere,  GITM  considers  the  densities  of  N2 
and  O2,  which  arc  the  main  constituents  at  100  km  altitude,  and  NO,  which  becomes  more  dominant  in  the 
upper  atmosphere. 

GITM  is  different  from  most  models  of  the  atmosphere  in  that  it  solves  the  full  vertical  momentum 
equation  instead  of  assuming  that  the  atmosphere  is  in  hydrostatic  equilibrium,  where  the  pressure  gradient 
is  balanced  by  gravity.  While  this  assumption  is  fine  for  the  majority  of  the  atmosphere,  in  the  auroral  zone, 
where  significant  energy  is  dumped  into  the  thermosphere  on  short  time-scales,  vertical  accelerations  often 
occur.  This  heating  causes  strong  vertical  winds  that  can  significantly  lift  the  atmosphere  [30], 

The  grid  structure  within  GITM  is  fully  parallel,  using  a  block-based,  two-dimensional  domain  decom¬ 
position  in  the  horizontal  coordinates  [70, 7 1],  The  number  of  latitude  and  longitude  blocks  can  be  specified 
at  runtime,  so  the  horizontal  resolution  can  easily  be  modified.  GITM  has  been  run  on  up  to  256  processors 
with  a  resolution  as  fine  as  0.31°  latitude  by  2.5°  longitude  over  the  entire  globe  with  50  vertical  levels, 
resulting  in  a  vertical  domain  from  100  km  to  roughly  600  km.  This  flexibility  will  allow  us  to  validate 
accuracy  by  running  data  assimilation  and  input  reconstruction  at  various  levels  of  resolution. 

It  is  possible  to  “fly”  satellites  through  GITM.  By  listing  the  times  and  locations  of  measurement  points, 
GITM  can  track  the  path  of  a  satellite,  outputting  simulated  data  at  the  specified  times  and  positions.  This 
feature  simplifies  implementation  and  validation  of  data  assimilation  and  input  estimation.  Fig.  5  shows 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


GITM  simulation  results  during  the  four  days  corresponding  to  the  CHAMP  data.  These  results  were  output 
every  minute  from  the  simulation  at  the  CHAMP  location  for  direct  comparison  with  the  satellite.  This  type 
of  comparison  allows  us  to  see  how  the  simulation  is  reproducing  storms  in  the  upper  atmosphere. 

Because  drivers  arc  important  to  the  upper  atmosphere,  GITM  can  be  driven  by  a  wide  variety  of  solar 
brightness  models  [45,74,95,100]  and  measurements  [20,98,99];  auroral  specifications  [40,42,69];  and 
high-latitude  electric  potential  models  [37,43,79,97],  In  addition,  the  assimilative  mapping  of  ionospheric 
electrodynamics  (AMIE)  technique  [13, 55, 56, 75, 76, 80]  can  be  used  to  more  accurately  estimate  the  high- 
latitude  aurora  and  electric  potential. 

5.2  Physics-Driven  Adaptive  Modeling  (PDAM) 

The  thermosphere  and  ionosphere  mark  the  true  transition  of  the  atmosphere  to  the  space  environment.  The 
top  of  the  thermosphere  (the  exosphere)  is  where  the  collision  frequency  of  the  neutral  particles  drops  low 
enough  that  it  becomes  difficult  to  describe  them  as  fluids.  In  the  lower  thermosphere,  the  neutrals  arc 
still  dominated  by  dynamics,  such  as  solar  tidal  forcing,  but,  as  altitude  increases,  and  the  density  of  ions 
also  increases,  the  neutrals  becomes  more  strongly  controlled  by  ion  forcing.  Also,  the  density  of  the  gas 
becomes  so  low  that  the  flow  speeds  can  become  quite  large,  reaching  over  1000  m/s  on  occasion.  This 
speed  is  comparable  to  the  sound  speed,  so  it  is  possible  to  get  supersonic  flows  in  the  upper  atmosphere. 
Furthermore,  because  the  flows  and  sound  speed  can  be  so  large,  the  dynamics  of  the  upper  atmosphere  arc 
truly  global — waves  propagate  from  one  side  of  the  Earth  to  the  other  in  just  a  few  hours.  This  means  that 
models  need  to  consider  the  global  system  and  not  just  focus  on  regional  dynamics. 

Ions,  because  they  arc  charged  particles,  arc  driven  by  electric  fields,  which  are  much  stronger  than  forces 
such  as  gravity  and  the  gradient  in  pressure.  The  electric  forces  are  so  strong  that  the  ions  quickly  reach  an 
equilibrium  velocity  that  can  have  a  significant  amount  of  small-scale  spatial  structure  [54].  The  neutrals, 
on  the  other  hand,  have  much  more  inertia  and  arc  slow  to  react,  so  the  ions  serve  to  transfer  momentum 
and  energy  (through  friction)  to  the  neutrals  [29,57],  This  transference  often  happens  on  small  scales  in 
regions  of  large  electron  densities  and  strong  electric  fields,  that  is,  in  the  auroral  zone  [31],  The  heating  on 
small  scales  causes  localized  perturbations  in  the  thermospheric  (neutral)  temperature  and  pressure,  which 
increase  winds  and  eventually  causes  a  global  disturbance  [12],  Therefore,  it  is  often  important  to  be  able  to 
capture  small-scale  dynamics.  The  problem  is  that  the  global-scale  phenomena  need  to  be  captured  as  well. 

Further  complicating  the  issue  is  the  fundamental  modeling  of  the  aurora  and  how  it  deposits  energy 
into  the  thermosphere  and  ionosphere.  The  aurora,  in  essence,  is  a  beam  of  electrons  that  is  shot  into 
the  thermosphere.  Many  researchers  have  simulated  how  the  atmosphere  reacts  to  this  beam  of  energy, 
and  parameterized  how  this  phenomenon  deposits  energy  into  the  system  [34,38].  However,  there  are 
aspects  of  this  process  that  are  often  simplified  and  arc  not  taken  into  account  when  simulating  active  time 
periods.  For  example,  models  of  the  auroral  energy  deposition  typically  use  a  static  thermosphere  that  is 
not  disturbed  [84],  The  composition  is  typically  static  in  the  models,  so  the  aurora  precipitates  into  the 
same  atmosphere  continually.  What  really  happens,  though,  is  that,  as  the  electrons  penetrate  into  the 
atmosphere,  heat  is  deposited,  preferentially  lifting  the  heavier  species  (such  as  N2  and  O2),  which  changes 
the  collision  frequency  between  the  beam  of  electrons  and  the  atmosphere  [17, 91, 92],  This  is  normally  not 
a  large  concern  since  the  amount  of  energy  deposited  is  not  substantial  enough  to  change  the  composition 
significantly,  but  during  large  auroral  events,  this  feedback  should  be  considered. 
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6  DDDAS  Frontier:  Mathematical  and  Statistical  Algorithms 

6.1  Context  and  Motivation  within  DDDAS 

The  ionosphere  and  thermosphere  arc  strongly  driven  by  the  Sun.  This  means  that  the  evolution  of  the 
thermosphere  depends  primarily  on  the  input  from  the  Sun  rather  than  its  current  state.  In  other  words,  the 
effect  of  an  initial  state  tends  to  be  “washed  away”  by  the  input.  Consequently,  estimates  of  solar  drivers 
can  greatly  enhance  the  accuracy  of  state  estimates  obtained  from  data  assimilation  with  unknown  inputs. 
Obtaining  accurate  estimates  of  these  drivers  represents  an  application  and  challenge  to  DDDAS  concepts 
and  technology. 

In  the  case  of  the  upper  atmosphere,  the  extreme  ultraviolet  radiation  produces  photo-ionization,  which 
in  turn,  through  chemistry  and  heating,  drives  the  properties  of  the  ionosphere  and  thermosphere.  Since 
a  significant  portion  of  the  EUV  and  X-ray  radiation  is  absorbed  in  the  atmosphere,  it  is  not  possible  to 
measure  the  flux  from  the  ground.  Instead,  a  proxy  is  used.  The  flux  solar  irradiance  at  a  wavelength  of  10.7 
cm  (called  F10.7)  is  thus  measured  on  the  ground  and  used  to  estimate  the  solar  spectrum  in  the  0-150  nm 
range.  The  unit  of  F  10.7  is  10~22  W/Hz/m2,  which  is  equivalent  to  1  solar  flux  unit. 

Although  good  estimates  of  F10.7  can  enhance  data  assimilation  in  the  ionosphere-thermosphere,  esti¬ 
mates  of  F10.7  arc  infrequent  (for  example,  once  per  day)  and  approximate.  This  situation  motivates  the 
need  to  estimate  F10.7  concurrently  with  the  states.  This  is  a  problem  of  input  estimation. 

6.2  Retrospective  Cost  Adaptive  Input  and  State  Estimation  (RCAISE)  for  DDDAS 

State  estimation  techniques,  such  as  the  Kalman  filter  and  its  numerous  variants,  use  measurements  to  re¬ 
cursively  refine  state  estimates.  In  the  simplest  case,  the  system  of  interest  has  the  form 

x(k  +  1)  =  Ax{k)  +  Bu{k)  +  w(k),  (6.1) 

y{k)  =  Cx(k)  +  v(k),  (6.2) 

where  x{k)  is  the  state  and  y(k )  is  the  available  measurement.  The  input  to  the  system  is  modeled  as  a 
combination  of  a  known  deterministic  signal  u(k)  and  an  unknown  stochastic  signal  w(k).  The  known 
deterministic  signal  u(k)  is  injected  numerically  into  the  observer  (not  shown  here),  which  uses  knowledge 
of  the  statistics  of  w(k)  and  the  sensor  noise  v(k)  as  well  as  knowledge  of  the  matrices  A  and  C  to  obtain 
an  estimate  x(k)  of  the  state  x  (/;:).  This  is  how  an  estimator  works.  In  practice,  however,  the  deterministic 
input  u(k)  is  often  unknown,  or  at  least  partially  uncertain.  A  common  practice  is  thus  to  treat  u(k)  as  paid 
of  the  stochastic  input  w(k).  This  approach,  however,  can  yield  poor  state  estimates  for  the  simple  reason 
that  the  characteristics  of  u(k)  may  be  quite  different  from  the  assumed  statistics  of  w(k). 

Because  of  uncertainty  in  u{k),  extensive  research  has  been  devoted  to  developing  extensions  of  the 
Kalman  filter  that  are  either  insensitive  to  knowledge  of  the  deterministic  input  or  that  attempt  to  estimate 
this  signal  in  addition  to  the  states.  These  techniques  arc  referred  to  as  unbiased  Kalman  filters,  unknown 
input  observers,  input  estimators,  and  state  estimators  with  input  reconstruction.  The  literature  is  extensive 
[14,36,49,59,60,72,94,96, 101],  but  no  unified  techniques  have  yet  emerged,  and  application  to  large- 
scale  data  assimilation  has  been  limited. 

The  importance  of  input  estimation  is  due  to  the  fact  that,  for  many  systems,  the  effect  of  the  initial  state 
decays,  and  thus  the  asymptotic  response  is  governed  entirely  by  the  forcing.  In  the  linear  case,  systems  with 
this  property  arc  asymptotically  stable.  In  the  nonlinear  case,  these  systems  arc  called  incrementally  stable 
or  contractive,  and  the  relevant  phenomenon  is  called  entrainment  [9,  81, 85].  The  presence  of  entrainment 
suggests  that,  at  least  in  some  cases,  the  accuracy  of  state  estimation  can  be  greatly  enhanced  by  the  ability 
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to  perform  input  estimation.  The  goal  of  this  research  on  data  assimilation  is  thus  to  reconstruct  the  inputs 
to  the  system  concurrently  with  the  estimates  of  the  state.  If  the  inputs  can  be  estimated  recursively,  then 
the  accuracy  of  the  input  reconstruction  can  be  significantly  enhanced  in  comparison  with  state  estimation 
alone. 

In  support  of  the  DDDAS  Frontier  Mathematical  and  Statistical  Algorithms,  this  task  is 
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Figure  5 :  The  atmospheric  density  p  at  407  km  altitude 
measured  by  the  CHAMP  satellite  (solid  line)  along  with 
simulation  results  from  GITM  (dashed  line).  The  up¬ 
per  plot  show  one-minute  values  of  the  results,  while  the 

lower  plot  shows  orbit-averaged  (i.e.,  90,-minute)  values.  .  .  . 

focusing  on  the  development  and  application  of  Retrospective  Cost  Adaptive  Input  and  State  Estima- 
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tion  (RCAISE)  [26].  The  RCAISE  architecture  is  shown  in  Figure  6,  where  the  error  z(k)  represents  the 
difference  between  measured  and  estimated  quantities  in  the  ionosphere-thermosphere.  RCAISE  is  an  out¬ 
growth  of  research  we  have  performed  on  adaptive  control  and  model  refinement  using  the  retrospective 


optimization  technique  described  in  [25, 28, 46-48, 65 

Within  the  context  of  data  assimilation  for  large- 
scale  applications,  input  estimation  can  be  performed 
by  applying  ensemble  Kalman  filter  techniques  [3- 
7,16,52,63].  These  techniques  employ  an  ensem¬ 
ble  of  simulations,  which  can  be  chosen  to  include  a 
distribution  of  representative  values  of  the  unknown 
input.  The  computational  burden  of  this  approach 
thus  increases  in  relation  to  the  number  of  ensem¬ 
ble  members.  Unlike  ensemble  filtering,  however, 
RCAISE  uses  a  single  simulation  model  as  the  ba¬ 
sis  of  both  state  estimation  and  input  reconstruction. 
Consequently,  unlike  particle  filters,  RCAISE  is  an 
ensemble-free  technique  for  data  assimilation  with 
input  reconstruction.  The  price  paid  for  this  conve¬ 
nience,  however,  is  the  lack  of  an  error  distribution 
on  the  estimated  input. 

The  prior  literature  on  state  estimation  with 
input  estimation  [14,36,49,60,72,94,96,101]  pro¬ 
vides  alternative  techniques  for  this  problem.  How¬ 
ever,  these  techniques  are  largely  confined  to  lin¬ 
ear  systems,  whereas  DDDAS  applications,  such 
as  space  weather  monitoring,  are  typically  highly 
nonlinear.  For  systems  with  nonlinear  dynamics, 
RCAISE  does  not  rely  on  an  analytical  model  of  the 
system,  requiring  only  the  ability  to  update  the  simu¬ 
lation  of  the  system,  and  this  simulation  can  be  in  the  : 


,82,93]. 


Figure  6:  State  and  input  estimation  architecture. 
Based  on  the  estimation  error  z  =  y  —  y,  the  adap¬ 
tive  algorithm  adjusts  the  the  feedback  block,  which 
estimates  the  unknown  input  u. 

orm  of  a  large-scale  computer  code,  such  as  GITM. 


6.3  Earlier  Research  Results 

We  now  summarize  recent  results  describing  the  application  of  RCAISE  to  the  ionosphere-thermosphere. 
An  overview  is  given  in  [24].  The  objective  of  input  estimation  supports  the  goals  of  DDDAS  by  enhancing 
prediction  accuracy.  In  particular,  state  and  input  estimation  are  considered  in  [2]  based  on  the  3D  GITM 
model  with  both  synthetic  and  real  measurements.  Here  we  consider  the  case  of  real  measurements  from 
the  CHAMP  satellite. 

This  study  is  motivated  by  the  fact  that  radio  propagation  and  satellite  drag  are  affected  by  the  Sun’s 
influence  on  the  ionosphere  and  thermosphere.  In  particular,  extreme  ultraviolet  (EUV)  and  X-ray  radiation 
produce  photo-ionization,  which,  in  turn,  through  chemistry  and  heating,  drives  the  formation  of  the  iono¬ 
sphere  and  shapes  the  thermosphere.  In  addition,  the  effect  of  the  EUV  and  X-ray  radiation  is  sufficient  to 
render  the  ionosphere-thermosphere  a  strongly  driven  system  [9,  81,  85]. 

Since  a  significant  portion  of  EUV  and  X-ray  radiation  is  absorbed  by  the  atmosphere,  it  is  not  possible 
to  measure  these  quantities  from  the  ground.  Instead,  a  proxy  is  used.  The  most  common  proxy  for  EUV 
and  X-ray  radiation  is  the  flux  solar  irradiance  at  a  wavelength  of  10.7  cm  (F10.7),  which  is  measured  (in 
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units  of  10“ 22  W  Hz  1  m“2  =  1  solar  flux  unit  (SFU))  by  the  Dominion  radio  observatory  in  Penticton, 
Canada  [68].  A  shortcoming  of  this  technique  is  that  F10.7  does  not  have  a  one-to-one  correlation  with  each 
of  the  wavelengths  in  the  EUV  and  X-ray  bands,  and  thus  the  measured  F10.7  is  often  a  misrepresentation  of 
the  true  solar  spectrum. 

Although  our  ultimate  goal  is  to  estimate  the  true  flux  in  multiple  EUV  and  X-ray  wavelength  bins,  a 
more  attainable  intermediate  goal  is  to  estimate  the  value  of  F10.7  that  best  characterizes  the  ionosphere  and 
thermosphere.  The  ability  to  estimate  F10.7  from  alternative  measurements  can  provide  a  cross  check  on  the 
available  measurements,  while  also  providing  an  illustrative  proof-of-concept  demonstration  of  RCAISE  as 
a  first  step  toward  estimating  X-ray  and  EUV  in  multiple  bands.  Furthermore,  current  models  do  not  fully 
capture  the  dynamics  of  the  ionosphere-thermosphere,  in  which  case  F10.7  can  be  used  as  an  input  to  the 
model  in  order  to  eliminate  the  errors  between  real  measurements  and  synthetic  (simulated)  measurements. 
This  study  thus  attempts  to  specify  F10.7  based  on  simulated  measurements  of  the  atmosphere  as  well  as 
with  real  satellite  data.  The  specified  F10.7  can  then  be  used  to  obtain  improved  estimates  of  the  state  of  the 
ionosphere  and  thermosphere  globally  and  possibly  predict  its  future  evolution.  This  is  a  problem  of  state 
and  input  estimation. 

To  estimate  F10.7,  we  use  the  Global  Ionosphere  Thermosphere  Model  (GITM)  [77].  GITM  simulates  the 
density,  temperature,  and  winds  in  the  thermosphere  and  ionosphere  across  the  globe  from  100  km  to  600 
km  altitude,  depending  on  the  solar  conditions  at  the  time.  The  main  inputs  to  GITM  are  the  high-latitude 
electrodynamics  (i.e.,  the  aurora  and  the  associated  electric  fields),  tides  from  the  lower  atmosphere,  and  the 
brightness  of  the  sun  at  various  wavelengths,  which  can  be  proxied  through  the  use  of  F10.7.  GITM  solves 
for  the  chemistry,  dynamics,  and  thermodynamics  of  the  upper  atmosphere  self-consistently  by  accounting 
for  interactions  among  various  species  of  ions  and  neutrals. 

In  [2]  we  used  the  retrospective  cost  adaptive  state  estimation  (RCAISE)  technique  [27]  to  estimate  the 
unknown  solar  driver  F10.7  using  the  Global  Ionosphere-Thermosphere  Model  and  satellite  measurements. 
RCAISE  assumes  that  the  input  to  the  system  is  unknown,  and  uses  retrospective  optimization  to  construct  an 
input  to  the  adaptive  estimator  that  minimizes  the  retrospective  cost  function.  The  retrospectively  optimized 
input  is  then  used  to  asymptotically  drive  the  error  between  the  measured  output  and  the  estimator  output 
to  zero.  In  this  way,  RCAISE  asymptotically  estimates  the  unknown  input  to  the  system  and  the  unknown 
states  of  the  system.  A  useful  feature  of  RCAISE  is  that  an  explicit  nonlinear  or  linearized  model  is  not 
required.  In  addition,  unlike  ensemble -based  data- assimilation  algorithms  [8,33,50],  RCAISE  uses  only 
one  copy  of  the  system  model  and  thus  is  ensemble-free. 

For  the  case  of  real  satellite  data  case,  the  neutral  mass  density  data  measured  by  CHAMP  (the  “truth 
data”)  is  labeled  y(k),  while  the  neutral  mass  density  data  measured  by  GRACE  is  labeled  yc(k).  First,  we 
run  GITM  with  the  measured  Fio.7(&)  and  record  the  neutral  mass  density  at  CHAMP  locations  and  GRACE 
locations,  which  are  labeled  ym(k)  and  yG,m{k),  respectively.  Next,  we  combine  RCAISE  and  GITM,  and 
use  y{k)  to  estimate  F 1  r j . 7 ( />' )  and  states.  The  neutral  mass  density  output  from  GITM  with  RCAISE  at 
CHAMP  and  GRACE  locations  are  labeled  y(k )  and  fjc(k),  respectively.  Note  that  data  from  GRACE 
are  used  only  as  a  metric  for  assessing  the  accuracy  of  state  estimates,  and  are  not  used  by  RCAISE.  We 
further  divide  this  setup  into  two  cases.  First,  in  RCAISE,  we  use  GITM  with  photoelectron  heating.  When 
photoelectron  heating  is  used  in  GITM,  then  the  neutral  density  output  from  GITM  at  CHAMP  locations 
using  measured  Fio.7(k’)  closely  matches  CHAMP  neutral  density  measurements.  However,  it  should  be 
noted  that  the  photoelectron  heating  efficiency  coefficient  that  yields  low  error  between  GITM  and  CHAMP 
is  obtained  by  trial  and  error,  and  cannot  be  calculated  or  measured.  In  the  second  case,  in  RCAISE,  we 
use  GITM  without  photoelectron  heating.  In  this  case,  GITM  with  measured  F10.7 (k)  yields  a  large  error 
between  the  outputs  from  GITM  at  CHAMP  locations  and  CHAMP  measurements.  In  this  case,  RCAISE 
will  use  Fi0.7(fc)  as  an  input  to  GITM  in  order  to  correct  the  errors  between  CHAMP  measurements  and 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


the  output  from  GITM  at  CHAMP  locations,  and  thus  account  for  the  inaccuracies  incurred  by  removing 
photoelectron  heating  from  GITM. 


(a)  Output-matching  performance  of  RCAISE  for  GITM  us-  (b)  State-estimation  performance  of  RCAISE  for  GITM. 
ing  driver  estimates. 

Figure  7:  (a)  /igo  ,y(k),  /j,go^(k),  pg9  ,ym(k),  &9o,y(k),  <190, y(k),  and  0-99, ym{k)  for  the  case  of  real  CHAMP  satellite 
data  and  GITM  with  photoelectron  heating.  For  this  example,  GITM  with  RCAISE  yields  6%  lower  RMS  (z)  compared 
to  GITM  with  measured  Fio.7(fc).  (b)  This  plot  shows  /j,90^yG(k),  /j,90}ya(k),  M90 ,yG,m{k),  090 ,ya{k),  <J9o,yG(k),  and 
a90,yG  m(fc)  for  real  GRACE  satellite  data  and  the  case  of  real  CHAMP  satellite  data  and  GITM  with  photoelectron 
heating.  For  this  example,  GITM  with  RCAISE  yields  11%  reduction  in  RMS(zg)  compared  to  GITM  with  measured 
F10.7  (k). 

Let  p(k)  G  M  be  an  arbitrary  signal,  and  let  T  be  a  positive  integer.  Then,  for  all  k  >  T,  define  the 
windowed  average  of  the  signal  p(k)  as 


VT,P{k)  =  ^  P(*)> 

i=k-T-\- 1 

where  T  is  the  interval  over  which  the  signal  is  averaged.  Similarly,  for  all  k  >  T,  define  the  windowed 
standard  deviation  of  the  signal  p(k)  as 


A 


VT,p{k)  = 


\ 


Y  (p(®)  -  ^T,p{i)y 

i=k-T-\-l 


The  root  mean  square  value  of  p(k)  is  denoted  by  RMS(p). 

When  GITM  is  used  with  RCAISE,  we  keep  Fio.r(fc)  at  a  constant  value  of  100  SFU  for  the  first  24  h, 
after  which  RCAISE  is  turned  on.  This  allows  the  response  due  to  initial  conditions  to  decay  significantly. 
We  implement  GITM  on  four  processors  with  a  resolution  of  5°  latitude  by  5°  longitude.  The  time  step  for 
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Figure  8:  Measured  and  estimated  Fio.7(A;)  for  the  case  of  real  CHAMP  satellite  data  and  GITM  with 
photoelectron  heating.  This  plot  shows  that  ylliQ  p  (A:)  (the  average  of  Fio.7(fc)  over  1  day)  converges  to 
within  6  SFU  (solar  flux  units)  of  the  measured  values  of  Fio.7(A;)  in  72  h. 

GITM  is  set  at  2  sec,  and  RCAISE  is  used  to  update  Fio.7(fc)  every  60  sec.  In  all  figures  in  this  section,  the 
vertical  black  line  indicates  when  RCAISE  is  switched  on.  Finally,  the  numerical  experiments  consider  the 
period  from  2002-11-24  to  2002-12-06. 

We  consider  the  case  where  the  truth  data  are  recorded  by  CHAMP  from  2002-11-24  to  2002-12-06  and 
we  use  GITM  with  photoelectron  heating.  We  calculate  the  RMS  of  z(k)  after  144  h  in  order  to  minimize  the 
effect  of  transients  in  y{k)  generated  during  the  convergence  of  the  adaptive  subsystem.  First,  we  consider 
the  measurements  from  CHAMP.  Figure  7(a)  shows  the  windowed  mean  and  variance  of  y(k),  y(k),  and 
ym(k).  For  this  example,  GITM  with  measured  Fio.7(fc)  yields  RMS(z)  =  6.5  x  10-13,  and  GITM  with 
RCAISE  yields  RMS  (2)  =  6.1  X  10-13.  In  other  words,  GITM  with  RCAISE  yields  6%  reduction  in 
RMS  (2)  compared  to  GITM  with  measured  Fio.7(&). 

Figure  8  shows  the  measured  and  estimated  Fio.7(A;).  This  plot  shows  that  ylliQ  p  (k)  (the  average  of 
Fi0.7(fc)  over  1  day)  converges  to  within  6  SFU  of  the  measured  values  of  Fio.7(fc)  in  72  h. 

Finally,  we  consider  data  from  GRACE  to  assess  the  quality  of  the  state  estimates.  Define  zc(k)  = 
Hoik)  —  yc{k)-  Figure  7(b)  shows  the  windowed  mean  and  variance  of  yc(k),  ycik ),  and  yG,m(k).  For 
this  example,  GITM  with  measured  Fio.7(fc)  yields  RMS(^g')  =  4  x  10-13,  whereas  GITM  with  RCAISE 
yields  RMS(zg)  =  3.6  x  10-13.  Therefore,  GITM  with  RCAISE  yields  11%  reduction  in  RMS(^g) 
compared  to  GITM  with  measured  Fio.7(&)- 


7  DDDAS  Frontier:  Application  Measurement  Systems  and  Methods 

7.1  Context  and  Motivation  within  DDDAS 

The  in-situ  measurements  required  for  GITM  are  difficult  to  make.  First,  they  must  be  made  in  space,  that  is, 
in  the  ionosphere  and  thermosphere.  Space  is  a  harsh  environment  where  extreme  thermal  fluctuations,  high 
vacuum  levels,  and  high-energy  radiation  require  specialized,  hardened  electronics  and  systems.  These  spe¬ 
cialized  space-based  sensors  are  expensive  to  build.  Also,  launches  are  infrequent  and  expensive.  Second, 
these  measurements  are  multiscale  in  time  and  space.  At  orbital  velocities  (~  8  km/s),  steep  gradients  in 
the  thermosphere  approach  rapidly  and  require  fast  sampling.  Our  measurement  environment  is  on  a  global 
scale,  covering  more  than  577  million  square  kilometers.  This  is  beyond  current  sensing  and  space  system 
capabilities. 
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A  potentially  desirable  feature  of  the  system  under  study,  the  thermosphere,  is  the  event-driven  nature 
of  the  system;  solar  storms  drive  thermospheric  changes.  Storm-time  dynamics  are  the  least  understood 
and,  consistent  with  the  goals  and  vision  of  DDDAS,  that  is  precisely  when  high  resolution  data  is  most 
valuable  to  the  data  assimilation  software.  Thus,  a  collaboration  between  the  applications  modeling  (GITM) 
and  the  system  operating  and  tasking  the  satellites  is  needed  to  task  the  sensing  systems  (satellites)  in  real 
time  in  order  to  adjust  sensing  capabilities  while  accounting  for  physical  constraints  on  the  satellites  such  as 
power  availability.  In  space  systems,  the  science  team  and  the  satellite  operations  team  arc  often  not  tightly 
coupled,  residing  in  separate  locations,  and  communication  between  them  is  slow.  It  can  take  days  to  weeks 
to  re-task  the  satellite  for  new,  advanced  operational  modes.  Thus,  our  effort  to  enable  real-time  tasking  of 
the  satellites  by  the  modeling  system  represents  a  novel  and  highly  nontrivial  advance  in  capability. 

The  foundation  for  our  space  weather  monitoring  effort  is  a  low-cost  satellite  constellation  composed  of 
Cubesats,  an  industry  standard  satellite  form  factor  that  has  a  space-qualified  launch  vehicle  interface,  the 
P-Pod  [15,67],  More  than  fifty  Cubesats  have  flown  on  rockets  from  the  US,  Europe,  and  Asia,  and  arc 
being  built  by  institutions  ranging  from  universities  to  Government  space  programs  (NASA)  to  established 
aerospace  corporations  such  as  Lockheed  Martin  and  Boeing.  At  the  University  of  Michigan,  we  have 
developed  and  launched  several  Cubesat  missions  for  space  weather  research  and  technology  demonstration. 
Developed  for  the  National  Science  Foundation,  two  satellites  for  the  Radio  Aurora  Explorer  (RAX)  mission 
were  launched  in  2010  and  2011,  respectively,  to  study  plasma  instabilities  that  lead  to  magnetic  field- 
aligned  irregularities  of  electron  density  in  the  lower  polar  thermosphere  (80-300  km)  [11],  RAX-2  is 
currently  fully  operational  and  supporting  bistatic  radar  experiments  with  transmitters  located  in  Alaska 
and  Canada.  We  are  also  developing  CADRE,  the  Cubesat-investigating  Atmospheric  Density  Response 
to  Extreme  driving,  to  make  measurements  specifically  for  thermospheric  modeling  and  data  assimilation. 
CADRE  will  be  ready  for  launch  in  late  2013  [22, 78], 

GITM  requires  neutral  densities  of  the  main  constituents  of  the  thermosphere,  in  particular-,  the  neutral 
winds,  the  ion  densities  of  the  main  ions,  and  the  ion  flows.  Instruments  are  available  for  obtaining  these 
measurements.  The  Ionospheric  Miniaturized  Electrostatic  Analyzer  (iMESA)  provides  electron  density  and 
temperature  measurements  using  a  miniaturized  electrostatic  spectrum  analyzer  [32],  The  Gated  Electro¬ 
static  Mass  Spectrometer  (GEMS)  is  based  on  an  electrostatic  energy  analyzer  and  provides  thermospheric 
composition  data  [44],  These  sensors  are  compatible  with  the  Cubesat  form  factor  [35], 

We  thus  have  a  challenging,  multiscale  sensing  requirement  from  our  advanced  model  system  that  is  not 
possible  with  current  technology.  The  foundation  for  this  sensing  system  is  a  constellation  of  approximately 
40  low-cost  Cubesat-based  sensors.  Cubesats  are  extremely  inexpensive  compared  to  traditional,  mono¬ 
lithic  satellites  and  thus  provide  an  enabling  technology  for  distributed  sensor  networks  in  space.  However, 
fundamentally  new  capabilities  are  needed  to  enable  this  system. 

7.2  Research  Effort:  Event-Based  Sensor  Reconfiguration  (EBSR) 

Within  the  DDDAS  Frontier  Application  Measurement  Systems  and  Methods,  we  are  developing  a  new 
capability  to  support  space-based,  multiscale  measurements,  namely,  event-based  sensor  reconfiguration 
(EBSR).  We  consider  sensing  platforms  and  support  systems  that  are  reconfigurable.  For  example,  dynamic 
elements  could  include  data-collections  rates,  sensor  selection,  sensor  physical  distribution,  and  sensor  di¬ 
rectionality.  Consistent  with  DDDAS,  sensor  reconfiguration  is  most  advantageously  based  on  events,  both 
physical  from  the  environment  and  feedback-based  from  the  data  assimilation  and  operations  system.  This 
will  enable  closed-loop,  real-time  sensor-system  modification  consistent  with  DDDAS  goals  and  vision.  Al¬ 
though  motivated  by  the  space-based  domain,  the  fundamental  characteristics  of  EBSR  can  be  summarized 
in  terms  of  three  primary  elements:  1)  the  system  must  be  able  to  reconfigure;  2)  the  system  must  be  able 
to  determine  when  to  reconfigure;  and  3)  the  system  must  collect  relevant  telemetry  and  data  to  inform  the 
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determination  system. 


Our  work  this  past  year  emphasized  in-field  sensor  calibration,  which 
addresses  elements  1)  and  3).  Calibration  is  the  development  of  a  transfer 
function  that  maps  sensor  readings  to  actual  (that  is,  physical)  values  as 
well  as  a  quantification  of  the  errors  in  the  sensor  measurements.  By 
definition,  reconfiguration  changes  systems  properties,  and  this  typically 
requires  sensor  recalibration.  Thus,  by  performing  in-field  calibration, 
we  can  extend  our  capabilities  for  reconfiguration. 

7.2.1  Magnetometer  Calibration 

We  finalized  our  method  for  in-field,  attitude-independent  magnetometer 
calibration,  which  includes  the  effect  of  time-varying  bias  due  to  elec¬ 
tronics  on-board  a  vehicle.  The  calibration  estimates  magnetometer  scale 
factors,  non-orthogonality,  and  both  constant  and  time-varying  bias.  The 
unique  contribution  of  this  effort  is  the  estimation  of  time-varying  mag¬ 
netometer  bias  due  to  nearby  electronics.  It  is  accomplished  by  includ¬ 
ing  vehicle  telemetry  in  the  measurement  model  and  estimating  constant 
parameters  that  map  time-varying  currents  to  magnetometer  bias.  Cali¬ 
bration  was  demonstrated  based  on  flight  data  from  the  RAX- 1  satellite, 
and  was  shown  to  significantly  reduce  the  uncertainty  of  off-the-shelf 
magnetometers  embedded  within  the  satellite  and  subject  to  spacecraft¬ 
generated  fields.  This  method  simplifies  the  satellite  design  process  by 
reducing  the  need  for  booms  and  expensive  magnetic  cleanliness,  result¬ 
ing  in  reduced  satellite  development  time  and  cost.  In  addition  to  on-orbit 
magnetometers,  the  calibration  method  is  applicable  to  air,  ground,  and 
water  navigation. 

The  developed  model  is  given  by 

r 

Bx  —  &BX  +  Xq  +  ^  ^  Sj,xlj  "F  Vxi  (7.1) 

2—1 

r 

By  =  b[By  cos (p)  +  Bx  sin(p)]  +  2/o  +  ^  +  Vy,  (7.2) 

2—1 


Figure  9:  The  Radio  Aurora  Ex¬ 
plorer  (RAX)  triple  Cubesat  devel¬ 
oped  at  the  University  of  Michigan. 


Bz  =  c[Bx  sin(A )+By  sin(^)  cos(A )+Bz  cos(fi)  cos(A)]+^o+X  ^zh+Vz- 


2—1 

(7.3) 

The  recalibration  procedure  estimates  magnetometer  scale  factors,  non¬ 
orthogonality,  and  constant  as  well  as  time-varying  bias.  The  time- 
varying  bias  is  captured  in  the  term,  )C[=]  Sijli,  j  6  {x.  y.  z},  where 
Sij  is  the  coefficient  that  maps  the  i-th  current  measurement  /,  to  the 

magnetic  field  in  the  j-th  magnetometer  axis,  and  r  is  the  number  of  current  measurements  included  in  the 
model.  Although  current  measurements  are  required  for  re-calibration,  on-board  current  sensors  are  typi¬ 
cally  part  of  spacecraft  health  monitoring.  Hence  inclusion  of  the  current  measurements  in  the  recalibration 
does  not  require  additional  sensors  for  the  purpose  of  recalibration.  Most  importantly,  this  model  does  not 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


require  any  knowledge  about  the  layout  of  the  electronic  components;  the  recalibration  requires  only  current 
measurements,  magnetometer  measurements,  and  the  expected  magnitude  of  the  geomagnetic  field. 

The  calibration  model  was  applied  to  data  from  the  RAX- 1  satellite  [23] .  From  the  data  presented  in 
Figure  7.2.1,  we  found  that  currents  produced  by  the  solar  panels  degrade  the  magnetometer  measurements. 
There  are  four  body-mounted  solar  panels  on  RAX-1,  and  we  include  measurements  of  the  current  in  each 
panel.  From  further  experimentation  with  the  on-orbit  data,  we  have  seen  that  a  fifth  element,  the  current 
drawn  from  the  battery,  degrades  the  measurements  [89].  Therefore,  we  include  these  five  current  measure¬ 
ments  in  the  magnetometer  model. 

The  results  of  the  calibration  arc  shown  in  Figure  1 1.  Figure  1 1(a)  shows  the  magnitude  of  the  corrected 
measurements  overlaid  with  the  expected  field  magnitude.  Figure  11(b)  shows  the  difference  between  the 
measured  and  expected  magnitudes,  and  Figure  11(c)  is  a  histogram  of  the  difference.  Figures  11(a)  and 
11(b)  show  the  significant  improvement  of  the  magnetometer  data  compared  to  the  corresponding  time- 
invariant  calibration  results  shown  in  Figures  10(b)  and  10(c).  The  root  mean  squared  error  (RMSE)  of  the 
measurements  after  calibration  for  time-invariant  errors  is  903  nT,  where  error  is  defined  as  the  difference 
between  the  expected  magnitude  and  the  magnitude  of  the  corrected  measurements.  After  calibration  with 
time-varying  bias,  the  RMSE  is  reduced  to  174  nT,  an  improvement  factor  of  5.2.  This  factor  corresponds 
to  an  order  of  magnitude  improvement  in  angular  accuracy. 

The  resolution  of  the  raw  magnetometer  measurements  is  128  nT  along  each  axis.  The  resolution  of  the 
corrected  sensor  readings  in  the  orthogonal  frame  is  transformed  by  the  scale  factors  and  non-orthogonality 
angles.  After  transformation  of  the  128  nT  resolution  with  the  calibration  parameters,  the  resolution  of  the 
calibrated  sensor  is  144,  143,  and  111  nT  along  the  x,  y,  and  z  axes,  respectively.  This  results  in  a  231 
nT  resolution  on  the  magnitude  of  the  corrected  measurements.  The  RMSE  of  the  calibrated  data  is  only 
174  nT,  which  is  below  the  sensor  resolution  of  221  nT  and  indicates  that  the  accuracy  of  the  calibrated 
measurements  has  approached  the  precision  limit  of  the  sensor. 

We  have  developed  a  method  for  on-orbit,  attitude-independent  magnetometer  calibration  that  includes 
time-varying  bias  due  to  nearby  electronics.  Existing  time-invariant  calibration  methods  have  been  extended 
by  including  current  measurements  in  the  sensor  model,  and  the  resulting  calibration  method  estimates  both 
time-invariant  errors  and  time-varying  bias  due  to  on-board  electronics.  The  time- varying  bias  is  estimated 
by  mapping  current  to  magnetometer  bias  through  constant  parameters.  The  calibration  parameters  are 
estimated  using  an  iterative  least-squares  minimization,  and  only  the  magnetometer  measurements,  current 
measurements,  and  the  expected  magnitude  of  the  geomagnetic  field  arc  required  for  the  calibration. 

The  usefulness  of  the  calibration  has  been  demonstrated  by  application  to  on-orbit  data.  In  application 
to  the  RAX- 1  satellite,  the  calibration  successfully  mitigates  magnetometer  bias  due  to  the  satellite  power 
system.  The  calibration  improved  the  RMSE  of  magnetometer  measurements  from  792  nT  to  219  nT, 
providing  an  order  of  magnitude  improvement  in  angular  attitude  accuracy. 

Traditionally,  time-varying  bias  is  mitigated  by  either  using  a  boom  to  physically  separate  the  magne¬ 
tometer  from  the  spacecraft,  or  by  using  design  and  manufacturing  practices  to  minimize  the  influence  of 
electronic  components  on  magnetometers.  Such  design  practices  increase  design  time  and  cost,  and  physical 
separation  of  a  magnetometer  from  on-board  electronics  may  not  be  possible  as  satellites  continue  to  de¬ 
crease  in  size.  The  algorithm  presented  in  this  work  effectively  replaces  such  design  practices  with  improved 
processing  of  the  sensor  measurements.  The  algorithm  utilizes  current  sensors  throughout  the  spacecraft, 
which  arc  sampled  at  the  same  time  as  the  magnetometers.  Therefore,  inclusion  of  current  sensors  and  the 
ability  to  sample  the  sensors  simultaneously  for  the  purpose  of  calibration  should  be  considered  in  the  de¬ 
sign  phase  of  the  vehicle.  Although  the  calibration  is  a  batch  method,  parameters  can  be  uploaded  to  the 
spacecraft  for  real-time  magnetometer  correction,  and  real-time  implementation  will  be  developed  in  future 
work.  The  calibration  has  been  applied  to  satellite-based  magnetometers  in  this  work,  but  this  algorithm  is 
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(a)  Magnitude  of  the  raw,  uncalibrated  measurements  (/rT)  overlaid  with  the 
expected  field  magnitude  using  the  IGRF  model. 


(b)  Magnitude  of  the  measurements  after  correcting  for  time-invariant  errors 
(pT).  The  differences  between  the  measured  and  IGRF  magnitudes  are  shown 
in  Figure  7.2.1. 


Minutes  Elapsed 

(c)  The  difference  (/rT)  versus  time.  The  sun  indicator  takes  the  value  of  one 
when  RAX-1  is  in  the  sun,  and  zero  when  in  eclipse,  which  shows  when  the 
solar  panels  are  illuminated  and  generating  current. 

Figure  10:  Data  from  the  RAX-1  PNI  MicroMag3  magnetometer.  The  .i-axis  of  each  plot  shows  time 
elapsed  since  the  start  of  the  data  set,  01-Dec-2010  08:30:46  UTC. 
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(a)  The  magnitude  of  the  corrected  measurements  (/rT)  versus  time  (minutes). 
The  magnitude  of  the  expected  magnetic  field  is  overlaid. 


(b)  Difference  between  the  corrected  measured  field  magnitude  and  the  ex¬ 
pected  magnitude  (/iT) 


(c)  Histogram  of  the  difference  between  the  measured 
and  expected  magnitudes  (nT).  There  are  5,405  total 
data  points. 

Figure  11:  Results  of  the  calibration  to  estimate  both  constant  errors  and  time-varying  magnetometer  bias. 
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applicable  to  other  magnetometer  applications  on  a  variety  of  platforms. 


7.2.2  On-Orbit  Calibration  of  Photodiodes  for  Attitude  Determination 

We  developed  a  method  for  on-orbit  calibration  of  photodiodes  for  sun  sensing  in  an  attitude  determination 
system  [90].  The  calibration  estimates  the  scale  factors  and  alignment  angles  of  the  photodiodes,  resulting 
in  higher  attitude  determination  accuracy  than  achieved  with  the  pre-flight  calibration  parameters.  The  cali¬ 
bration  is  implemented  with  an  extended  Kalman  filter  to  simultaneously  estimate  spacecraft  attitude  and  the 
calibration  parameters.  This  approach,  as  opposed  to  an  attitude-independent  method,  enables  the  calibra¬ 
tion  of  an  arbitrary  number  of  photodiodes  mounted  in  any  orientation  on  the  spacecraft  and  facilitates  the 
use  of  an  attitude-dependent  Earth  albedo  model.  The  method  is  demonstrated  by  application  to  flight  data 
from  the  RAX-2  satellite  and  results  in  an  average  angular  improvement  of  10°  in  sun  vector  measurements 
with  the  photodiodes.  Attitude  determination  accuracies  of  below  1°  in  each  axis  arc  demonstrated  using 
the  calibrated  photodiodes  in  combination  with  a  low-cost  three-axis  magnetometer  and  rate  gyroscope. 


8  DDDAS  Frontier:  Advances  in  Systems  Software 

8.1  Context  and  Motivation 

Previously,  we  discussed  the  need  for  new  methods  for  collecting,  analyzing,  and  using  data  gathered  from 
the  ionosphere-thermosphere  to  monitor  space  weather  and  thus  better  predict  the  orbital  motion  of  large 
space  objects.  To  successfully  meet  the  goals  outlined  in  these  three  sections,  we  also  need  to  develop  new 
models  for  scheduling  and  communicating  policies  on  how  to  transmit  these  data  to  Earth. 

Although  we  have  posed  (and  will  continue  to  discuss)  this  problem  in  the  context  of  Cubesat  data 
collection,  the  results  that  we  are  developing  have  much  broader  applicability.  Specifically,  we  consider 
problems  in  which  we  have  multiple  clients  (e.g.,  Cubesats)  and  multiple  servers  (e.g.,  ground  stations). 
The  clients  arc  largely  controlled  by  a  central  manager  who  provides  policies  to  guide  behavior;  the  clients 
in  turn  make  local  decisions  based  on  these  policies  as  well  as  on  environmental  factors  that  they  experience 
(for  example,  a  Cubesat  might  gather  more  data  than  originally  planned  if  significant  variation  between  its 
expected  and  observed  location  were  detected).  The  central  manager  has  discrete  and  possibly  infrequent 
opportunities  to  communicate  with  the  clients  to  update  the  policies.  Finally,  the  clients  arc  dependent  on 
the  servers,  which  themselves  may  be  controlled  by  external  agents,  may  have  conflicts  with  other  users, 
and  may  be  subject  to  uncertainty  in  their  availability  and  their  performance  characteristics.  Problems  of 
this  type  can  be  found  in  computing  systems,  telecommunication  systems,  energy  systems,  UAV  systems, 
and  many  more  real-world  contexts.  To  address  these  fundamental,  challenging,  and  widely  applicable 
scheduling  problems,  we  arc  developing  heuristics,  optimization-based  algorithms,  and  simulation  tools  for 
hierarchical  closed-loop  scheduling  (HCLS).  Our  work  towards  these  goals  is  described  below. 

8.2  Modeling  Space  Communication  Networks 

We  have  developed  a  general,  analytical  framework  for  modeling  an  operational  satellite  mission.  We  define 
a  framework  as  a  set  of  reusable  elements  and  templates  for  describing  dynamics,  constraints,  and  goals. 
The  framework  analytically  represents  the  dynamic  interaction  of  states  (such  as  position,  energy,  and  data) 
and  subsystem  operations  (such  as  communication  and  energy  management)  of  an  operational  satellite.  It 
captures  mission  constraints,  which  arc  often  called  requirements,  that  specify  minimum  performance  levels. 
It  also  enables  analytical  expression  of  objectives,  which  arc  goals  of  the  mission  to  be  max-  imized.  The 
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framework  is  generic  and  modular  such  that  it  is  capable  of  supporting  a  variety  of  mission  architectures  and 
scenarios. 

8.2.1  Model  Elements 

The  four  main  elements  of  the  framework  arc  parameters,  states,  subsystems,  and  the  schedule.  Elements 
arc  constant  or  time-dependent,  time  notation  is  omitted  for  simplicity. 

Parameters  -  A  parameter,  p,  is  a  model  input  that  provides  numerical  values  to  dynamically  model  sys¬ 
tem  states  and  subsystem  functions.  Let  P  be  the  set  of  all  model  parameters,  where  p  6  P.  Examples 
parameters  arc  orbital  parameters,  ground  station  locations,  and  Tf. 

States  -  A  system  state  is  a  model  variable,  and  is  defined  as  the  information  at  some  initial  time  that, 
combined  with  the  input  (parameters  and  the  schedule)  for  all  future  time,  uniquely  determines  the  output  for 
all  future  time  [21].  Let  X  =  [x\, ....  a;/-,  ...xm}T  be  the  vector  of  all  the  system  state  variables,  where  there 
are  m  variables.  Example  states  include  on-hoard  resources  such  as  energy  and  payload  data.  Opportunities 
for  mission  operations  such  as  payload  operation  and  ground  station  availability  arc  also  system  states. 
An  opportunity  is  modeled  as  binary,  o  G  {0, 1},  where  a  value  of  one  indicates  an  opportunity  and  zero 
indicates  no  opportunity. 

Subsystems  -  A  subsystem,  s,  performs  functions  on  states.  Let  S  be  the  set  of  all  subsystems.  A  single 
function  operates  on  state  k  and  is  denoted  fs,j,k  £  A,  where  j  6  Js  is  the  function  index,  and  Js  is  the  set 
of  all  function  indices.  fs,j,k  is  an  element  of  the  set  of  functions,  F. 

Schedule  -  The  schedule,  U ( t ),  is  a  series  of  time-dependent  events  that  describes  how  and  when  the  sub¬ 
system  functions  operate  on  the  states.  Events  are  scheduled  when  there  arc  opportunities.  For  example,  a 
data  download  event  may  occurs  when  there  is  a  line  of  sight  between  a  ground  station  and  satellite.  The 
schedule  is  designed  to  achieve  the  mission  objectives  while  satisfying  the  mission  constraints.  The  schedule 
may  be  an  output  (e.g.  when  a  solver  is  used  to  find  an  optimal  schedule)  or  an  input  (e.g.  when  simulating 
a  given  schedule  to  test  performance). 

8.2.2  Framework  Formulation 

The  model  is  formulated  as  a  conventional  optimization  problem  in  Eqs.  8. 1-8.4.  Mission  objectives, 
represented  in  Eq.  8.1,  maximize  the  total  transfer  of  a  mission-specific  system  state,  x*,  a  component  of  X, 
over  the  planning  horizon.  The  decisions  in  the  optimization  problem  arc  when  and  how  the  events  occur, 
which  arc  captured  in  the  schedule,  U (t),  which  is  an  output  of  the  optimization  problem  as  formulated 
here.  The  constraints  in  the  formulation  include  state  dynamics  (Eq.  8.2),  bounds  on  state  values  (Eq.  8.3), 
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and  mission  requirements  (Eq.  8.4). 


ma x{x*(Tf)  } 

U(t )  J 

s.t. 

X(t  +  At)  =  N(X,  P,  t)  +  ^2  E  F*j(X’  u>  *) 

s£S  j^:Js 


X min  —  X(t)  —  Xmaa; 
Ai+i 

0k, i  <  /  x:fe(t)  dt 
Jti 


(8.1) 

0  <  t  <  Tf  (8.2) 

0  <  t  <  Tf  (8.3) 

Vxfc  G  X,  i  G  Ik  (8.4) 


States  evolve  over  time  due  to  nominal  dynamics  and  subsystem  functions  (see  Eq.  8.2).  Nominal 
dynamics  arc  independent  of  subsystem  functions.  The  vector  of  nominal  dynamics  equations  is  defined  in 
Eq.  8.5,  where  each  element  k  represents  the  nominal  dynamics  of  state  Xk-  Orbital  motion  and  battery  self¬ 
discharge  arc  example  nominal  dynamics  of  the  state  variables  position  and  on-board  energy,  respectively. 

N(X,P,t)  =  [ni(X,P,t),...,nfc(X,P,t),...]T,  (8.5) 


The  vector  of  subsystem  functions  that  operates  on  the  state  vector  is  expressed  in  Eq.  8.6.  The  inputs  to 
each  function  fs,j,k  include  the  states,  parameters,  schedule,  and  time.  Note  the  vector  in  Eq.  8.6  contains 
zero  entries  when  combined  subsystems  and  functions  do  not  operate  on  specific  states. 


F aj(X,U,PaJ,t)  =  [fsjA(X,U,P8j,t),.„,fsjtk(X,U,PaJ,t),...]T  Vs  G  S,j  G  Js 


(8.6) 


The  nominal  and  functional  dynamics  in  Eqs.  8.5  and  8.6  may  each  be  described  by  any  type  of  function, 
for  example  they  may  be  analytical  or  extracted  from  a  simulation  system. 

The  state  vector,  X,  is  constrained  by  lower  and  upper  bounds,  {Xmjn,  Xmax}  G  P,  as  in  Eq.  8.3. 
Example  bounds  include  maximum  and  minimum  battery  capacity  and  maximum  data  storage  capacity. 

Operational  mission  requirements  arc  represented  in  Eq.  8.4  by  enforcing  a  minimum  change  in  system 
state  over  a  specific  time  period.  For  example,  there  may  be  a  mission  requirement  that  a  minimum  amount 
of  state  (such  as  energy)  must  be  acquired  or  consumed  during  a  certain  period  of  time.  Each  interval  i  G  /),., 
where  Ik  is  the  set  of  intervals  spanning  the  full  planning  horizon  for  state  X}-,  has  a  start  time,  0  <U<Tf, 
where  the  end  of  interval  i  corresponds  to  the  start  of  interval  i  +  1.  Eq.  8.4  enforces  a  minimum  change 
of  state  Xk  during  every  interval  i  G  Ik,  represented  as  0k, i-  The  change  in  state  during  interval  i  is  its 
integrated  time  rate  of  change  from  t,  to  tl+\.  For  states  without  requirements,  0k, i  will  be  zero  Vi  G  Ik- 

Another  perspective  for  describing  spacecraft  operations  is  to  consider  subsystem  functions  individually. 
In  particular,  consider  the  analytical  relationship  between  inputs  and  outputs  specific  to  subsystem  s  and 
function  j,  Zsj  =  gs,j{ YSJ-,  U.  P,  t),  where  the  vector  of  inputs  is  Ysj  and  the  vector  of  outputs  is  Z.SJ, 
which  are  both  comprised  of  components  of  X.  The  function  gsj  is  the  combination  of  fs,j,k  VA:  G  K ,  i.e., 
it  models  the  impact  of  subsystem  s  and  function  j  on  all  state  inputs  and  outputs.  These  relationships  for  a 
single  subsystem  arc  shown  in  Figure  12. 


8.2.3  Block  Diagram  Representation 

We  represent  the  model  framework  using  a  conventional  control  system  block  diagram  to  demonstrate  the 
interaction  of  the  various  model  elements  in  Figure  13.  The  set  P  of  parameters  is  provided  to  the  input 
block,  which  identifies  opportunities  for  subsystem  functions,  O,  and  interprets  the  mission  requirements, 
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Figure  12:  A  generic  representation  of  the  subsystem  function  Zsj  =  gsj(Ysj,U,  Psj,t)  for  s  =  1  and 
j  =  1,  2,  3.  All  values  are  time  dependent. 


R,  as  control  inputs.  The  error  signal  is  expressed  as  E  =  R  M.  where  M  is  estimated  state  values,  which 
are  measured  by  on-board  or  ground  sensors.  E,  P  and  R  are  provided  to  the  scheduler,  which  generates 
the  operational  schedule,  U.  Note  that  U  is  an  output  of  the  controller  and  an  input  to  the  dynamic  system. 
The  states  evolve  according  to  both  the  nominal  dynamics  and  subsystem  functions  as  prescribed  by  the  U, 
where  updated  states  (after  time  At)  are  denoted  X(t  +  At).  Unmodeled  realistic  disturbances,  D,  may 
be  injected  into  the  system  and  modify  the  state.  Mission  performance  is  evaluated  by  measuring  the  states 
and  verifying  if  the  mission  requirements  are  satisfied  and  comparing  realized  objectives  to  their  expected 
values.  Feedback  control  occurs  when  the  scheduler  updates  U  according  to  mission  performance,  i.e.,  uses 
E  in  for  future  scheduling  decisions. 
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Figure  13:  Elements  and  dynamics  of  the  system  model  represented  with  a  conventional  feedback  control 
loop  diagram.  The  non-italicized  labels  are  the  conventional  elements  of  a  control  feedback  loop.  The 
italicized  labels  are  the  elements  of  the  modeling  framework. 


8.3  Research  Plans 

Thus  far  we  developed  a  sophisticated  model  for  space  communication  systems  and  developed  an  opti¬ 
mization  algorithm  for  single-satellite,  multiple  ground  station  scheduling.  Our  next  steps  will  included  the 
following: 
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Develop,  model,  implement  and  analyze  a  stochastic  version  of  the  single  satellite,  multiple  ground 
station  problem.  Initial  results  will  be  reported  later  in  the  year. 

Develop  concepts  for  multi-satellite  and  multi-ground  station  architectures. 
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